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FOREWORD 


The  multivariate  splines  presented  in  this  report  have  been  designed 
for  the  direct  analysis,  smooth  representation,  and  efficient  storage  of 
data  which  is  randomly  spaced  in  multiple  dimensions  and  which,  additionally, 
may  be  periodic  along  one  or  more  of  these  directions.  Of  special  interest 
to  the  U.S.  Naval  Oceanographic  Office  is  the  applicability  of  these  splines 
to  randomly  sampled  geophysical  and  oceanographic  features  and/or  processes. 

The  univariate,  minimum-curvature  spline  has  been  applied  as  a  smoother 
and  "gap  filler"  in  the  SEASAT  ephemeris  correction;  this  spline  also  has 
been  used  to  generate  a  pseudomonthly  Generalized  Digital  Environmental  Model 
(GDEM)  and  employed  to  model  the  north  wall  of  the  Gulf  Stream. 

Possible  future  applications  of  the  multivariate,  minimum-curvature 
splines  are  as  a  closed-form  representation  of  the  gravimetric  geoid  as  an 
efficient,  closed-form  upward  continuation  of  gravity  anomalies  and  as  a 
parametric,  closed-form  representation  of  GDEM  based  on  latitude,  longitude, 
depth,  and  time  of  month. 

This  work  was  carried  out  within  the  Advanced  Technology  Staff  and  has 
been  reviewed  by  its  head.  Dr.  Thomas  M.  Davis. 


Captain,  *SN 
Acting  Commanding  Officer 
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ABSTRACT 


This  paper  develops  a  collection  of  multivariate  splines  which  are  appli¬ 
cable  to  randomly  spaced  data.  The  splines  are  constructed  from  a  linear  com¬ 
bination  of  spline  basis  functions  having  evenly  spaced  nodes.  A  pseudoinverse 
solution  of  coefficients  is  obtained  via  a  singular  value  decomposition  (SVD) . 
Spurious  oscillations  are  damped  by  zeroing  out  an  appropriate  number  of  the 
singular  values.  A  transformation  is  performed  prior  to  entering  the  SVD  so 
that*  whenever  the  solution  is  nonunique,  the  "minimum-curvature"  solution  is 
obtained.  An  added  feature  of  these  splines  is  an  option  through  which  a  periodic 
fit  may  be  produced  in  one  or  more  dimensions  if  it  is  so  desired. 
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INTRODUCTION 


Geophysicists  and  oceanographers  are  constantly  faced  with  the  task  of 
analyzing  and  interpreting  randomly  spaced  data.  Moreover,  it  is  not  unusual 
for  such  data  to  have  been  collected  nonun if ormly;  some  portions  of  the  sampled 
area  may  be  data  dense,  while  other  portions  may  contain  no  data  at  all.  This 
nature  of  geophysical  sampling  presents  a  problem  to  the  analyst  attempting  to 
represent  the  data  as  a  parameterized  curve  or  surface.  What  often  happens  is 
that  the  linear  system  resulting  from  the  parameterization  has  an  infinite 
number  of  solutions,  or  is  highly  ill-conditioned,  or  both. 

In  this  report,  we  will  develop  a  collection  of  multivariate  splines 
which  are  immune  to  these  difficulties.  This  will  be  accomplished  by  first 
representing  the  splines  as  a  linear  combination  of  spline  basis  functions. 

The  linear  system  resulting  from  this  parameterization  will  then  be  transformed 
into  a  new  reference  frame  and  a  singular  value  decomposition  in  conjunction 
with  pseudoinversion  will  be  used  to  solve  the  new  system.  The  reference  frame 
transformation  will  be  derived  so  that  whenever  the  solution  is  nonunique,  the 
"minimum-curvature"  solution  is  obtained. 

Another  important  aspect  of  geophysical  and  oceanographic  data  is  that  of 
period '.city.  This  is  the  case  when  one  considers  geophysical  phenomena  on  a 
global  scale,  or  the  seasonal  dependence  of  many  types  of  oceanographic  data, 
the  latter  data  being  at  least  quasi-periodic  as  a  function  of  time  of  year. 

For  this  reason,  we  will  also  consider  the  concept  of  a  periodic  spline  and 
derive  an  elegant  way  of  representing  multivariate,  periodic  splines. 

The  first  step  in  the  development  is  the  construction  of  sets  of  basis 
functions  for  the  multivariate  splines. 


2.  MULTIVARIATE  SPLINE  BASIS  FUNCTIONS 

Let  S  be  the  set  of  all  cubic,  univariate  spline  functions  having  evenly 
spaced  nodes  (or  knots)  s^,  82*  •  •  • »  sn  separated  by  the  distance  h.  A 
function  f  is  a  member  of  S  if  an  only  if  f  is  a  cubic  polynomial  on  each  sub¬ 
interval  J"Sj ,  sj+]^  an<*  twice  continuously  differentiable  throughout  the 
interval  £s^  s^  . 

It  is  well  known  that  the  set  S  constitutes  an  (n+2) -dimensional  linear 
space  of  functions.  Thus  any  set  of  n+2  linearly  independent  members  of  S,  e.g. 
{  <f>  ;  j  =  1,  .  .  .  ,  n+2  }  ,  forms  a  basis  for  S.  This  means  that  each  f  £  S 
can  be  uniquely  represented  as 

n+2 

fit)  =  ^  *  ^(o.  (i) 

j-i 

One  of  the  most  commonly  used  bases  [Lawson  (1974),  Schoenberg  (1973), 
Ahlberg  (1967),  deBoor  (1978)]  can  be  defined  using  the  two  cubic  polynomials 

Pj  (0)  =  • 2503  (2) 


p2  (0)  -  1-. 75 (1+9) (1-0)  . 


With  0j  defined  as 


where 


J*1 ,  .  .  . ,  n— 1 


h  =  s  j  ,  J=l»  •  •  •»  n— 1 


the  following  set  of  functions, 


for  s^  <_  t  <  s._3  and  5  <  j  £  n+2 


P1  ^j-3^  for  sj-3  —  1  —  sj-2  and  ^  —  J  —  n+2 

P2  for  sj-2  —  t  —  sj-l  and  3  £  J  £  n+1 

P2  (l-0j_p  for  sj_^  £  t  <  Sj  and  2  <  j  <  n 

(I-®,)  for  s.  <  t  <_  s.  .  and  1  <  j  £  n-1 

0  for  s...  <  t  <  s  and  1  <  j  <  n-2 

1+1  —  —  n  —  J  — 


j  =  l . n+2 
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represents  a  basis  for  S.  This  spline  basis  is  illustrated  in  figure  1  for 
the  case  h=l,  n-11.  (Note:  Figures  are  located  on  pages  39  through  59.) 


Multivariate  spline  bases  can  be  generated  from  this  univariate  spline 
basis  through  the  use  of  tensor  products  [de  Boor  (1978)  and  Schultz  (1973)  J  . 
To  simplify  the  multivariate  development,  it  will  be  required  that  the  nodal 
spacing  in  each  direction  be  equal  to  the  distance  h.  This  can  be  done  with¬ 
out  any  loss  of  generality;  different  nodal  spacings  in  each  of  the  multi¬ 
directions  can  be  achieved  by  appropriate  scaling  of  the  independent  variables 
of  the  data. 

The  bivariate  spline  basis  is  composed  of  the  functions 


4>ij (ti «  t2)  =  <j)i(t1)^(t2) 


i=l ,  .  . 
j=l,  •  • 


n.+2 

(7) 

n2+2 


Orthogonal  subsets  of  this  basis  are  depicted  in  figures  2-17  for  the  case 
h=l,  n^=ll,  n2=ll.  The  entire  basis  is  represented  by  the  superposition  of 
these  figures. 

The  trivariate  spline  basis  is  made  up  of  the  functions 

i*l,  .  .  .,  nj+2 

*ijk(tl»  t2’  t3)  =  <^i(tl)<^j  (t2)<^k(t:3)  j=1’  *  ’  •»  V2  (8) 

k=l,  .  .  .,  n3+2 

In  general,  the  M-variate  spline  basis  consists  of  the  functions 

jj=l,...,n  j+2 

t2,”,,tM)  =  <i>j1(tl)<i!j2(t2)-”<f)jM(tM)  J2=1 . n2+2  (9) 

« 

JM=1, . . .  ,1^+2 


Hence,  any  bivariate  spline  with  evenly  spaced  nodes  can  be  uniquely  re¬ 
presented  as 


f( tj,  t2) 


i=l  j-1 


Xij^i(tl)<{lj(t:2)’ 


(1C 


any  trivariate  spline  with  evenly  spaced  nodes  can  be  uniquely  represented  as 


f(tj»  t2,  t3) 


n.+2  n_+2  n«+2 

tit- 


ijk<f>i(tl^j(t2)<!>k(t3)  ’ 


(11) 


and  so  forth. 
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3. 


THE  LINEAR  SYSTEM 


Using  such  bases,  the  problem  of  finding  a  multivariate  spline  that  best 
fits  a  set  of  data  reduces  to  solving  a  linear  system 


Ax  =  b 


(12) 


in  some  sense  (usually  a  least-squares  sense)  where  A  is  an  m  x  N  matrix.  For 
example,  in  the  univariate  case,  given  a  set  of  data  {(t^jy^):  t^e  [s^>sn]> 
i=l,...,m|  ,  employment  of  equation  (1)  leads  to  the  linear  system 

n+2 


3=1 

n+2 

j=i 


(13) 


n+2 

y  =  >  x  <J>  (t  ) . 

m  /  .J  j  j  m7 

3  =  1 


In  this  instance,  N  =  n+2  and 


w 

4>2(t1) 

'  '  •  <*n+2(tl) 

"X1 

’  yl  " 

A  = 

4>1(t2) 

4>2(t2) 

'  *  •  <|)n+2(t2) 

t  x= 

X2 

• 

,  b= 

y2 

• 

• 

*> 

♦iV 

w 

'  •  •  <J>n+2(tm) 

• 

_Xn+2  _ 

• 

(H) 


In  the  multivariate  cases,  use  of  equations  (10)  and  (11)  leads  to  similar 
linear  systems.  The  matrix  A  will  always  contain  the  basis  functions  evaluated 
at  the  appropriate  points  in  time  or  space,  the  vector  x  will  contain  the  spline 
coefficients  to  be  solved  for,  and  the  vector  b  will  contain  the  dependent  values 
of  the  data. 

If  the  data  is  relatively  uniformly  dense  and  there  exist  several  data 
points  within  each  nodal  division  (subinterval,  subarea,  subvolume,  etc.),  then 
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the  columns  of  the  matrix  A  will  be  linearly  independent  and  the  system  will 
have  a  unique  least-squares  solution 

=  A#  b*  (15) 

where 

A  =  (AXA)  A1  (16) 

T  -1 

is  the  pseudoinverse  of  A.  In  practice,  (A  A)  is  rarely  computed,  but  rather 
^  T 

x  is  obtained  by  performing  a  Cholesky  decomposition  on  A  A  or  employing  a  se¬ 
quential  Householder  processing  procedure  [Lawson  and  Hanson  (1974)]  . 

If  the  data  is  not  uniformly  dense,  and  there  exist  nodal  divisions  con¬ 
taining  no  data,  then  the  columns  of  the  matrix  A  may  be  linearly  dependent, 

T 

thus  making  A  A  singular.  Hence,  there  will  be  an  infinite  number  of  least- 
squares  solutions.  The  system  may  be  highly  ill-conditioned  as  well.  Under 
these  conditions,  the  Cholesky  and  Householder  approaches,  per  se,  break  down, 
and  an  alternative  approach  must  be  taken.  Such  an  alternative  is  the  use  of 
a  singular  value  decomposition  in  conjunction  with  pseudoinversion  (LeSchack 
(1976),  Golub  (1965),  (1970),  Van  Loan  (1976),  Rao  (1971),  Ben-Israel  (1974)]. 
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4.  THE  SINGULAR  VALUE  DECOMPOSITION 

The  singular  value  decomposition  of  an  m  x  N  matrix  A  of  real  elements 
[LeSchack  (1976) ,  Golub  (1965),  (1970),  Van  Loan  (1976)]  is  a  factorization 
of  the  form 

A  =  USkVT  (17) 

where  U  is  an  m  x  m  orthogonal  matrix,  V  is  an  N  x  N  orthogonal  matrix,  and  S 

K 

is  an  m  x  N  matrix  whose  only  nonzero  elements  are  K  positive  diagonal  elements 

h 

X^  ,  i  =  1,,  .  .,  K  called  the  singular  values  of  A.  The  singular  values  are 
arranged  in  descending  order  of  magnitude,  with  K  being  the  rank  of  A. 

Remarkably  efficient  algorithms  exist  to  numerically  perform  the  decompo¬ 
sition  (17).  An  SVD  algorithm  in  conjunction  with  pseudoinversion  represents 
a  practical  and  accurate  technique  for  solving  ill-conditioned,  least-squares 
problems  which  originate  in,  or  are  transformed  to,  the  form  of  equation  (12). 
Using  a  matrix  outer  product  expansion,  equation  (17)  can  be  written  as 


K 


xi  Vi’ 


(18) 


i-1 


,th 


where  u  and  v  are  the  i— =  column  vectors  of  the  matrices  U  and  V,  respectively 
h  1  th 

and  X^  is  the  i=  singular  value  of  the  matrix  A.  The  pseudoinverse  of  A  is 


then  given  by 


K 


J  m  ST*  ,-Js  -  T. 

A  "  Xi  "i  "i  * 


(19) 


i-1 


In  the  event  that  K=N,  this  pseudoinverse  is  identical  to  the  pseudo inverse 
(16).  For  K<N,  the  matrix  A  is  rank  defective,  and  the  pseudoinverse  (19) 
yields  the  least-squares  solution  of  minimum  length. 

Ill-conditioning  of  the  system  originates  from  some  singular  values  being 
relatively  much  smaller  than  other  ones.  This  ill-conditioning  can  be  overcome 
by  setting  the  smaller  singular  values  to  zero,  that  is,  terminating  the  summa  ■ 
tion  in  (19)  short  of  K.  Hence,  a  truncated  pseudoinverse  can  be  obtained  via 

<  ■  ]C  ^  ’A1  <2° 

i-1 


9 


where  L<K.  Whenever  L<K,  the  full  solution  (L=K)  has  been  deemed  unstable 
and  has  been  replaced  by  a  stable  solution  at  the  expense  of  a  slightly  larger 
least-squares  residual  error. 

A  procedure  for  determining  the  termination  index  L  was  presented  by 
Shim  and  Cho  (1981)  for  use  in  image  restoration  problems.  An  optimal  termina¬ 
tion  index  was  derived  assuming  that  both  the  signal  and  noise  were  white  (with 
respect  to  a  Karhunen-Loeve  transform)  and  is  given  by 

L  =  max  |k  :  X.  >  —  - —  \  (21 

opt  k  (  k  -  N-  (SNR)  ) 

where  SNR  is  the  signal-to-noise  ratio  and  X^  is  the  k=  singular  value  squar¬ 
ed.  This  simple  formula  produced  good  image  restorations  in  Shim  and  Cho's 
computer  simulations.  In  practice,  this  author  has  found  that  the  formula 
works  equally  well  in  curve-fitting  applications. 


The  SVD  technique  just  discussed  yields  the  least-squares  solution  of 
minimum  length.  Thus,  if  the  method  were  applied  directly  to  the  linear  system 
(12),  the  vector  of  spline  coefficients  would  be  minimized.  Since  it  is 
known  [deBoor  (1978)]  that  these  coefficients  model  the  surface  they  rep¬ 
resent  (i.e.,  a  plot  of  the  coefficients  yields  a  shape  similar  to  the  sur¬ 
face,  at  least  to  within  a  scale  factor),  minimizing  the  spline  coefficients 
is  roughly  equivalent  to  minimizing  the  spline  surface  itself.  Unfortunately, 
there  is  no  physical  justification  for  doing  so.  It  is  more  desirable,  there¬ 
fore,  to  produce  a  solution  which  optimizes  some  physically  justifiable  per¬ 
formance  index.  This  can  be  achieved  by  representing  the  system  in  a  new  ref¬ 
erence  frame  prior  to  application  of  the  SVD. 

Let  the  performance  index  P  be  given  by 

P  =  xTQx‘  (22 

where  Q  is  a  positive  definite  matrix  and  let 

T_* 

y  -  Lx  (23 

where  L  is  the  lower  triangular  Cholesky  factor  of  Q.  This  implies  that 


.  -1nT*. 

(L  )  y- 


(24 


Substituting  (24)  into  equation  (12),  we  have 

By^  =  b 

where 


(25a 


-1  T 

B  =  A(L  )  . 


(25b 


If  the  SVD  technique  is  now  applied  to  the  system  (25),  the  vector  y  will  be 
minimized,  forcing  the  performance  index 


_  vr  _  T-».  ,tT*,T/tT-  -JU 

P  =  x  Qx  =  x  LL  x  =  (L  x)  (L  x)  =  y  y 


(26; 


to  be  a  minimum.  The  appropriate  spline  coefficients  can  then  be  obtained  via 
equation  (24). 


The  problem  now  is  to  find  a  positive  definite  matrix  Q  such  that  P  will 
represent  a  physically  important  aspect  of  the  solution.  Before  doing  so,  how¬ 
ever,  the  concept  of  a  periodic  spline  is  introduced. 


6.  PERIODIC  SPLINES 

For  a  univariate  spline 


n+2 

5 


to  be  periodic  over  the  interval  £ a snj,  the  spline  along  with  its  first  and 
second  derivatives  must  match  at  the  interval  end  points,  that  is. 


f(s i>  “  f(sn) 

f'(s i>  -  r<sn> 


f"( sj)  -  f^(sn). 


This  means  that 


n+2 

£ 

n+2 

xjV8i>  *  £ 

Xj*j<8n) 

j-l 

j-i 

n+2 

n+2 

£ 

*j*i  ^i’  ■  £ 

*j*j'  <8n^ 

j-l 

j-i 

n+2 

n+2 

£ 

j-l 

vr  <8i>  ■  £ 

j-i 

Vj"  (8n> 

These  equations  reduce  to 

3 


£  IjV8l)  ■  £  Vj<8>>) 

j-l  j-n 

3  n+2 

<si:>  ■  £  'v 

j-n 

3  n+2 

£  vr  '“i*  ■  £  vr  (sn> 

j-l  j-n 


PREVIOUS  PAGE 
IS  BLANK 


since  these  are  the  only  nonzero  contributors  to  the  sums. 


Evaluating  the  functions  at  the  Interval  end  points  yields 


1  _i_  .  1 

4  Xi  +  X2  +  4  X3 

1  .1 

4  n  n+1  4  n+2 

(31a) 

"3  .  3 

4  XI  +  4  x3 

-3  3 

"  4  *1,  +  4  xn+2 

(31b) 

3  ,.3 

2  xi  -  3x2  +  j  x3 

‘  T  xn  -  3xn+l  +  f  xn+2. 

(31c) 

If  we 

let 

a 

‘  xl-xn 

(32a) 

b 

“  x2"xn+l 

(32b) 

c 

“  x3“Xn+2 

(32c) 

then 

system  (31)  becomes 

1 

4  a 

+  b  +  j  c  -  0 

(33a) 

-3 

4  8 

+  |c  -  0 

(33b) 

3 

2  a 

-  3b  +  |  c  -  0 

(33c) 

which 

has 

the  unique  solution  a*0. 

b*0,  c“0.  Hence,  the  equations 

X1  "  Xn 

(34) 

x2  =  Xn+2 

(35) 

x3  “  xn+2 

(36) 

represent  the  necessary  and  sufficient  conditions  for  f (t)  to  be  periodic  over 
the  Interval  [Sj,  snJ  [deBoor  (1978)]  . 

One  way  to  Implement  these  conditions  would  be  to  auguent  system  (12)  with 
the  three  additional  equations  (34)  -  (36) ,  strongly  weighting  these  latter  equa¬ 
tions;  that  would  only  slightly  increase  the  size  of  the  system.  However,  our 
aim  is  to  develop  a  collection  of  multivariate  splines,  and  in  the  higher 
dimensions,  the  number  of  conditions  which  must  be  satisfied  for  periodicity 
significantly  increases.  This  is  due  to  the  fact  that  the  multivariate  spline, 
along  with  its  first  and  second  partial  derivatives,  must  match  everywhere 
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along  the  boundaries  of  the  domain  (bivariate  area,  trivariate  volume,  etc.) 
considered.  Hence,  such  an  augmenting  procedure  is  rejected. 


Instead,  we  will  develop  an  elegant  and  efficient  way  of  representing 
and  computing  the  periodic  splines.  This  will  be  done  by  determining  a  spline 
basis  for  the  univariate,  periodic  spline.  Multivariate,  periodic  spline  bases 
can  then  be  generated  through  the  use  of  tensor  products.  By  taking  this 
approach,  the  size  of  the  system  to  be  solved  will  be  reduced. 


Substituting  equations  (34)  -  (36)  into  equation 

(27)  yields,  for 

n>4 , 

n-1 

3-1 

(37) 

where 

^(t)  -  4>j(t)  +  4>j+n_i(t) 

j=l ,  •  •  • ,  3 

(38a) 

^j(t)  =  ^(t) 

j“4 . n-1 

(38b) 

For  n-3 

xl  -  x3  =  x5 

(39a) 

x2  -  x4 

(39b) 

so  that 

f  (t)  »  Xjif^t)  +  2 ( t) 

(40) 

where 

^j(t)  -  ^(t)  +  <f>3(t)  +  cf>4(t) 

(41a) 

i^2(t)  *>  4>2(t)  +  cf>4 (t) . 

(41b) 

For  n«2 


so  chat 


f(t)  -  (t) 

where 

^(0  -  (^(t)  +  <}>2(t)  +  <j>3(t)  +  <|>A(t). 


(43! 


(44; 


Combining  these  results,  we  have  that 

n-1 

/(t)  x^(t) 

>1 


(45] 


where 


n+2 


V‘> 


i-j 
in  steps 
of  n-1 


(t) 


j-1 . n-1  (46) 


Hence,  the  set  of  functions  (t) ,  j»l,.  .  n-1 1  represents  a  basis  for  the 

(n-1) -dimensional  space  of  all  cubic,  univariate,  periodic  splines  with  evenly 
spaced  nodes.  Figure  18  depicts  the  first  three  basis  functions  ^(t),  ^(t) , 
i|>3(t)  for  the  case  h-1,  n-11.  (Note  their  periodicity.)  Figure  19  shows  a 
typical  univariate,  periodic  spline  fit. 

A  basis  for  a  bivariate  spline  that  is  periodic  in  both  the  tj  and  t2 
directions  is  given  by 

pu(v  v  . v; 

j“l»  •  •  •»  n2~ 1 


A  basis 
is  given  by 


for  a  trivariate  spline  that  is  periodic  in  all  three  directions 


i* 1,  .  .  .,  n^-1 

Pijk(tr  V  t3)  "  ^iCtl)^j(t2)\(t3)  *  *  *»  n2_1(48) 

kal,  .  .  .,  n^-I 


If  we  are  only  concerned  with  periodicity  in  one  direction,  e.g.,  the  t2 
direction,  then  appropriate  bicubic  and  tricub ic  spline  bases  can  be  generated 


T 


T 


T 


-k",T!L  ""J.  "I  ■  «  "  U  ' 


and 


Pij(tr  t2)  =  <j>i(tl)^j(t2) 


Pijk(tl’ 


’ 2 * 


t3) 


wwvv 


i-1, . . . »n  j+2 

j"  1 » •  •  •  *n2~ 1 

i=l . n^+2 

j=l » • • • »n2~l 
k*  I)***}  2 


(49) 


(50) 


Thus,  using  this  technique,  multivariate  splines  which  are  periodic  along 
one  or  more  arbitrary  directions  can  be  produced. 


The  elegant  part  of  this  construction  is  that  the  periodic  framework  now 
fits  neatly  within  the  regular  spline  framework.  Hence,  the  periodic  option 
will  require  no  additional  computer  storage.  Furthermore,  the  same  computer 
routine  that  evaluates  the  regular  spline  can  be  used  to  evaluate  the  periodic 
versions.  This  is  made  possible  by  re-representing  the  periodic  spline  in 
terms  of  the  regular  spline  basis.  For  example,  in  the  univariate  case,  once 
the  n-1  coefficients  have  been  obtained,  we  can  set 


Xn  =  X1 

(51a) 

xn+l  =  X2 

(51b) 

xn+2  "  x3 

(51c) 

then 

once 

represent  the  periodic  spline  via 
the  (n^-1)  x  (n2~l)  coefficients 

equation 

have  been 

(27) .  In  the 

obtained,  we 

bivariate 

can  set 

case. 

Xi,n 

"  xi  i 
2  1,1 

(52a) 

K 

Xi,n2+1 

"  Xi,2 

i-1,  .  . 

• ,  n^-1 

(52b) 

• 

and 

Xi,n2+2 

Xi»3 

(52c) 

X  4 

n  ^ ,  j 

"  Xi>j 

(53a) 

Xnj+1 , j 

X2,j 

j“l»  .  .  . 

,  n2+2 

(53b) 

\+2,i 

X3,j 

(53c) 

-a 
.  * 

.  * 
—I 

*! 


: 

■j 

•j 


j 


uO 

ti 
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then  represent  the  periodic  spline  with  equation  (10) .  In  the  trivariate 
case,  once  the  (n^-1)  x  (^-1)  x  (n^-l)  coefficients  have  been  obtained,  we 
can  set 


X  .  = 

i, j .n3 

Xi,j,l 

(54a; 

i, j .n^+l 

Xi,j,2 

i=l ,  .  . 

•  9  n 1 

(54b 

i,j,n3+2 

Xi,j,3 

j-1,  •  • 

•  9  iu-l 

(54c 

Xi,n2,k 

Xi,l,k 

(55a 

Xi,n3+l,k 

Xi,2,k 

i=l ,  .  . 

•  9 

(55b 

'i,n3+2,k 

Xi,3,k 

k*l ,  .  . 

n3+2 

(55c 

Xn, ,j,k 

Xl,j,k 

(56a 

tij+1,  j  ,k  =  X2,j,k 


xn1+2,j,k  "  X3,j,k 


j=l,  •  •  •  ,  n^+2 


k=  1 ,  .  •  . ,  n^ 2 


then  represent  the  periodic  spline  through  equation  (11).  If  periodicity  had 
been  enforced  in  only  one  direction,  e.g.  the  t^  direction,  in  the  two  previous 
examples,  then  steps  (52)  and  (53)  would  be  replaced  with  only  step  (53)  and 
steps  (54) -(56)  would  be  replaced  with  only  step  (56). 


7. 


"CURVATURE"  MATRICES 


One  physically  important  aspect  of  the  solution  is  that  of  curvature. 
Hence,  we  shall  proceed  to  develop  a  performance  index  that  is  a  measure  of 
the  "curvature"  of  the  spline  surface.  We  use  quotes  here  because  what  we 
will  actually  do  is  construct  a  performance  index  which  measures  the  Frobenius 
norm  of  the  second  derivative  of  the  spline  over  the  domain  of  interest.  This 
approximately  measures  the  curvature  over  that  domain  for  small  curvatures. 

In  general,  we  will  consider  quantities  of  the  form 

/h4  M| | V2f ( t1 ,  t2>.  .  .,tM) | |2dt1dt2.. .dt^  (5 
*'aM 

2 

where  V  f(t^,  t2>  ...»  t^)  is  the  second  derivative  (second  gradient)  of  the 
spline,  ||  •  ||  is  the  Frobenius  norm,  and  M  is  the  number  of  variates. 

The  form  (57)  reduces  to 

J  *  "x^Cx  (5i 

where  x  is  the  vector  of  spline  coefficients  and  C  is  a  positive  semidefinite, 
"curvature"  matrix.  The  factor  h"'  where  h  is  the  node  spacing,  is  simply  a 
scale  factor  which  has  no  effect  on  the  minimization  of  J.  However,  it  does 
serve  an  important  purpose;  its  presence  results  in  the  curvature  matrix  C 
being  a  function  of  the  grid  size  only,  and  not  of  the  grid  spacing  as  well. 
Since  C  is  only  positive  semidefinite,  an  epsilon  technique  will  be  used  to 
produce  the  positive  definite  matrix  Q  required  for  the  transformation  (23). 
The  matrix  Q  will  be  given  by 


Q  -  C  +  El  (51 

where  I  is  the  identity  matrix  and  £  is  a  small  number,  large  enough  only  for 
the  computer  at  hand  to  identify  Q  as  being  positive  definite.  Thus  the  per¬ 
formance  index  P  will  be  a  measure  of  the  weighted  sum  of  the  "curvature"  and 
the  length  of  the  vector  of  spline  coefficients.  Since  e  will  be  extremely 
small,  essentially  only  the  curvature  will  be  measured. 


7.1  UNIVARIATE  CURVATURE  MATRICES 


In  the  univariate,  nonperiodic  case, 


11  t 

fw  -  E)  Vj' 


The  second  derivative  of  f( t)  is  given  by 


\  •.  •«  \  \ 


so  that 


d2f(t) 

dt2 


E  vr<o 

j-i 


Jrb=s  r 

r  “  ■ 

a=Sj  L 


dt 

dt 


(62a) 


J^b  n+2  2 

1,3  E  vr(t>  dt 

a  U"1 


rb  ["  n+2  n+z 

■I  »3  E  vr<*>  E  vr(t>  dt 

,/a  i-l  J  LJ-l 


(62b) 


(62c) 


n+2  n+2 


i-l  j-1 


x±  h  3J  <p''(t)<p''(t)dt  Xj 


(6  2d) 


-jkT 

■  x  Wx 


(62e) 


« V.  *  -V  ‘  •>- 
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i 


where 


W(n) 


( Vn)) 


i= 1 ,  .  .  . ,  n+2 
j=l,  .  .  . ,  n+2 


(63) 


with 


Wij(n)  =  ^  fa  ^>i'( t)^(t)  dt. 


(64) 


The  weighting  matrices  W(n) ,  n=2,  3,  .  .  .  were  computed  by  integrating 
(64)  analytically;  n  was  increased  until  a  pattern  was  recognized  enabling 
generation  of  W(n)  for  arbitrary  n.  These  weighting  matrices  are  given  in 
appendix  A. 

In  the  univariate,  periodic  case, 


fit) 


n-1 

j-i 


V'>- 


The  second  derivative  of  f(t)  is  given  by 


so  that 


d2f(t) 

dt2 


n-1 

Lvr 

j=i 


(t) 


(65) 


•-1 

‘.’i 


(66)  ^ 


I 

v 

■j 

•j 
■ j 

jj 

8 


4 


j  =  r=s n  h3(" 

•4«s, 


1  2 


dt 


dt 


/b  rn-1 

1,3  Z/Zj" 

j“l 


-.2 


(t) 


dt 


(67a) 


(67b) 
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f 


r  fc,fg  *i*h  \t 

J< a  L^1  J  U"1 


(t)  dt 


(67c) 


=  ICS  xi  kf  *r<t>«r(t)dtlxj 

i-i  j=i  l  a  j 


(67d) 


where 


WP(n) 


("« (n) ) 


i=i ,  .  .  . ,  n-1 
k= 1$  •  •  • f  n- L 


h3  £  <^(t)^(t)dt. 


(67e) 


The  superscript  P  denotes  periodicity. 

The  weighting  matrices  W1* (n)  ,  n=2,  3,  .  .  .  were  computed  as  in  the  non¬ 
periodic  case,  integrating  (69)  analytically  and  establishing  a  pattern  for  the 
generation  of  (n)  for  arbitrary  n.  These  weighting  matrices  are  given  in 
appendix  B 

In  both  of  the  above  cases,  the  curvature  matrix  C  is  identical  to  the 
weighting  matrix: 


for  a  nonperiodic  spline  and 


for  a  periodic  spline. 


C(n)  =  W(n) 


C(n)  =  (n) 


Figure  20  shows  an  application  of  the  univariate,  miniraum-crvature 
spline  which  has  resulted  in  an  improved  NAVOCEANO  product.  A  gross  estimate 
of  the  north  wall  of  the  Gulf  Stream  is  obtained  from  a  gridded,  partially- 
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+h  J  ♦i1(tl)*j.(tl)dtl hJ  ♦i,(t2)*j"(t2)dt2  xj.j. 

•'a^  1  J1  •'a  2  2  J2  J1  . 


(74b) 


iij+2  11^+ 2  n^+2  n^+2 


V1  V1  V*  V1 


t±1i2  [Vi1J1(ni)Zl232<n2)  + 


+  2i131("l>“i232(n2>]  x3l: 


(74c) 


where  (n)  is  defined  as  in  (64)  and 


^(n)-h  f  <p'(t)<p'(t)dt 


zlj(n)m:tX  *±(t)*j(t)dt 


The  weighting  matrices 


(yu<n)) 

(■«w) 


i=l . n+2 


j“li  •  •  •»  n+ 2 


i“l . n+2 


J“1 . n+2 


i“l,  .  .  n+2 


j»l . n+2 


i=l ,  ,  n+2 


j“l . n+2 


were  computed  by  integrating  (75)  and  (76)  analytically  and  are  given  in 
appendices  C  and  D,  respectively. 

Setting 

%l  <V  n2>  '  Wi1j1(nl)zi232(n2)  +  2)'i1J1(nl)yi2j2(n2> 

+  zi1J1("l,wl2J2(n2>  <79> 


ft 

I 


I>  V  •  . 

vWU.Vj.1; 


k  -  +  (i2~l) (nj+2) 


•  •  > 


and 


l  -  jj  +  (J2-l) (nj+2) 


1 


V1, 


1 

n2+2 


j ^*1 »  •  •  • »  n 2 
j  2= 1 »  •  •  • »  n2+ 2 


(80a) 


(80b) 


N  =  (ni+2)(n2+2) 


(81) 


we  have 


H  N 

J  ‘  2  1’  ^ 
k-1  l  *1 

v>p  ^ 

*  X  C(n1>  n2)x. 


(82a) 


(82b) 


Here,  the  curvature  matrix  C(n^,  n2)  Is  defined  via  (79)  and  the  vector  x' 
contains  the  spline  coefficients  stored  columnwise. 


Figure  21  displays  an  application  of  the  bivariate,  minimum-curvature 
spline.  This  spline  was  applied  to  nineteen  data  points  identifying  Eddy  F 
in  the  Gulf  of  Mexico.  The  figure  depicts  the  resulting  minimum-curvature 
surface  and  its  contour. 


In  the  bivariate,  periodic  case. 


uj  i  2 

f(tl’  t2)  "  2  S  Xj1j2^1(tl)^j2(t2)’ 

V1  V1 


The  second  partials  of  f( t^,  t2)  are  given  by 


(t,) 


(84a) 


32f 

3t,9t 


lut2 


hi2>i 


u,)r'(t) 
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so  that 


dt^dt 


2 


P  P  p  p  “I 

+  2y .  .  (n,)y .  ,  (n„)  +  z,  ,  (n.)w.  .  (n0)  x,  , 


i-1 


>  • 


n-1 


where  w 


ij 


(n)  is  defined  as  in  (69)  and 


(n) 


^(tWj(t)dt 


J-l. 


n-1 


(n) 


^±(t)^  (t)dt 


i-1. 


J-l. 


.  ,  n-1 

. ,  n-1 


The  weighting  matrices 


YP(n)  =  (yij(n)) 

and 

ZP(n)  -  ^zij  <n^ 


i=l . n-1 


j-l . n-1 


i-1 . n-1 


j“l . n-1 


were  computed  by  integrating  (86)  and  (87)  analytically  and  are  given  in 
appendices  E  and  F,  respectively. 

Setting 


Ck,l  (V  V 


Wi  j  (nl)zi  j  (n2} 
1J1  2J2 


where 


k  “  ii  +  (v^^r1) 


V1* 


V1* 


• »  n^-1 

(91a) 

*  *  n2“^ 


j >  •  • 


•  9  n^-1 


*.  * 
** 


*  -  J,  +  (j 9~1) (n.-l) 


(92) 


N  =  (nj-1) (n2~l) 


we  have 


2  )xJ 


k-1  i  -1 


x  C(n^,  n2)x 


with  the  curvature  matrix  C^,  n2)  being  defined  by  (90) 


If  periodicity  is  desired  only  in  the  t^  direction,  then 


n^-1  n2+2 


*  ■  2Lf  f*  \h\ 
V1  V1 


<tlH32<t2) 


and  the  curvature  matrix  C  is  obtained  from 


Ck,!(V  n2}  "  Wi1j1(nl)zi2J2(n2)  +  2yi1j1(nl)yi2J2(n2) 


(93a) 


(93b) 


(94) 


+  z  (n  )w  (n  )  (95) 

±ljl  1  V2  2 


where 


i^=l ,  •  •  •  »  n^  1 


k  -  i1  +  (i2-l)(nrl) 


(96a) 


i2~  1 ,  •  •  §  *  n2"^ 2 


. nl"*^ 


l  -  j,  +  (J,-l)(n.-l) 


(96b) 


If  periodicity  is  desired  only  in  the  t2  direction,  then 


f( tj,  t2) 


n^+2  n2~l 


2-,  L,  X31J2*31(Il)*32(t2) 

3J-1  ~ 


and  the  curvature  matrix  C  is  obtained  from 


‘Tc.l  (ni*n2)  -  wi131(nl)*i2j2<n2)  +  2yi131(nl)l'l232(n2>  +  zi131<nl)wl232<” 


where 


k  -  i1  +  (i2-l) (nL+2) 


l  -  Jj  +  (J  2“ 1  >  (n1+2) 


ij  =  l »  •  •  •,  n 2 

i2=l»  •  •  • ,  n2~l 
j ,  n^+2 


j2*'1  * 


n2~l 


7.3  TRIVARIATE  CURVATURE  MATRICES 
In  the  trivariate  case, 


f (ti»  t2,  t3)  - 


"jlj2j35jl(tl>5j2(t2)5j3(t3) 


where 


5.  (t.)  •  4> .  (t.)  for  nonperiodicity  in  the  t.  direction 
Ji  1  Ji  1  x 

5,  (t.)  »  ip.  (t.)  for  perodicity  in  the  t.  direction 
Ji  Ji  1  1 

and 

d^  »  n^+2  for  nonperiodicity  in  the  t^  direction 
d^  ■  n^-1  for  periodicity  in  the  t^  direction. 


(97) 

2>  (98) 

(99a) 

(99b) 

(100) 

(101a) 

(101b) 

(102a) 

(102b) 
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Taking  the  second  partials  of  7(t^,  t2>  t^)  and  proceeding  as  in  the  two 
previous  sections. 


nr 


a3  a2  al 


/  2  \  2  /  2  \  2  /  2  \  2 

h  (<£_)  +  (it)  + 

>ti  /  \*i  /  VS  / 


/  A  A  +  2  &-Y 

VVW  V8tl3t3/ 


,9t2St3 


2_,  2^XkCk,l(nl’  V  nj 


f<=l  l  =1 


dtldt2dt3 


(103a) 


(103b) 


x  C(n1>  n2>  n3)x‘ 


where  the  curvature  matrix  C(n1>  n2>  n3)  is  given  by 


12  3 

K  ,  N  IT  ,  .  KJ 


(103c) 


12  3 

IT  \  KT  .  K"3 


C?c,l  <nl’  n2’  n3)  =  Wi1j1(nl)zi2j2(n2)zi3j3(n3)  +  zi1j1(nl)wi2j2(n2)zi3j3(n3) 


12  3 

K1  s.K 


12  3 

K1  %  K  ,  jr 


+  zi1j1(nl)zi2J2(n2)wi3j3(n3)  +  2yi1j1(nl)yi2j2(n2)zi3J3(n3) 

12  3 

+  2y^  .  (n.)z^  (n,)y^  .  (n3) 

iljl  1  *232  2  i3J3  3 


+  2ziljl(nl)yi2j2(n2)yi3j3(n3) 


(104) 


K1  =  "blank"  for  nonperiodicity  in  the  t^  direction 
Kd  =  "P"  for  periodicity  in  the  t^  direction 


(105a) 

(105b) 


i^l,  .  .  dl 


k  =  i.  +  (i  -l)d  +  (i  -l)d.d  i2-l,  .  •  d  (i06a) 


(106b) 


I  -  Jl  +  (j2"1)d1  +  (j3-l)d1d2  j2-l.  .  .  d2 
and  j  ^=1 )  i  •  *  t  ^ 

N  =  djd^.  (107) 

P  P  P 

The  weighting  matrices  w^  (n)  ,  (n)  ,  (n)  ,  w^  (n)  ,  (n)  ,  z  (n)  are 

defined  as  in  (64),  (77),  (78),  (69),  (86),  (87),  respectively. 

For  example,  if  periodicity  is  desired  only  in  the  direction,  equa¬ 
tion  (104)  becomes 


CM(nl-  V  V  ■  ”l1j1(nl>2i2j2<n2)2i3J3<"3)  +  2iljl(nl)wi2j2(n2)2i3j3(n3) 


+  2l1J1(,1l)2i2j2(n2)»i3J3("3)  +  2>'i131<nl)’'i232(n2)2i3j3(n3) 


+  2>'i1j1<nl)2l232(n2)>'i3j3<n3) 


+  22i131(nl)5i232(:’2)yl333C"3)- 


(108) 


7.4  QUADRIVaRIATK  curvature  matrices 

In  the  cuadrivariate  case. 


f  (ti»  ^2*  t3*  = 


(109) 


with  E,  (t.)  and  d  ,  i«l, 2,3,4  defined  as  in  (101)  and  (102).  Taking  the  second 
Ji  1  1 


partials  of  f( t^>  t2»  t^,  t^)  and  proceeding  as  before, 


'k4  [l3  rb 2  fb  1 


J  - 


-'a^  s a  */a2  •'a^ 


m  -m -m 


a2f  X2 


3t, 


+  2 


(kk) 


2 


2 


(-&-Y+  2(-&-) 

W'J  \9t29t3/ 


2  +  2  (^J  +  2  (n$J  «iW*4 


(110a) 


N  N 


*kck,l  (nl’  n2’  n3*  n4)xI 


k-1  1=1 


(110b) 


-  x  C(n1,  n2,  n3,  n^)x 


(110c) 


where  the  curvature  matrix  C(n^,  n^.  n^,  n^)  is  given  by 


%,l  (V  n2*  “S’  V  ’  "i1J1(l,l,2i2i2<n2)zl3j3<n3>2i4J4("4> 

+  2<n2>  %  j3("3)  %J4  <V 

12  3  4 

+zi131("l>Zi232(n2)wi3J3<n3)2l4J4<”4) 

+  zi,  3  J  (nl>  zi2j2("2>  zl333'”3>\34  (n4> 

12  3  4 

+2yl131<n1>yi2j2<n2)zi333("3>zi4J4(”4) 


+  2yijJ  J  <nl>zI23  2<n2>yl333  (n3>  %J4  <n4> 

+2yi131(nl)zi232(n2)zl333(n3)yi4j4(n4> 

+  2zi131(al>yi232("2)yl333("3)zi434(n4) 

+  2zl131(nl>yi232<n2>zl333("3)yi434<n4) 
+  2^131<nl)zy2<"2)yi333<"3>yi434(n4) 


(111) 
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8.  SUMMARY  AND  CONCLUDING  REMARKS 


A  collection  of  multivariate,  minimum-curvature  splines  for  randomly 
spaced  data  has  been  developed.  The  splines  are  represented  as  a  linear  com¬ 
bination  of  spline  basis  functions  having  evenly  spaced  nodes.  This  para¬ 
meterization  results  in  a  linear  system  of  the  form 

Ax  =  "b  (115) 


where  the  m  x  N  matrix  A  contains  the  spline  basis  functions  evaluated  at  the 
appropriate  points  in  time  or  space,  the  vector  x  contains  the  spline  coeffi¬ 
cients  to  be  solved  for,  and  the  vector  b  contains  the  dependent  values  of  the 
data.  A  reference  frame  transformation  is  then  performed  to  yield 


A(L  “  b. 


(116) 


The  matrix  L  is  the  lower  triangular  Cholesky  factor  of  a  positive  definite 
matrix  Q  which  essentially  is  a  measure  of  the  curvature  of  the  spline  surface. 
The  matrix  Q  is  given  by 


Q  *  C  +  El  (117) 

where  C  is  the  so-called  "curvature"  matrix,  1  is  the  identity  matrix,  and 
E  is  a  small  number,  large  enough  only  for  the  computer  at  hand  to  identify 
Q  as  being  positive  definite.  The  new  system  (116)  is  then  solved  using  a 
singular  value  decomposition  in  conjunction  with  pseudoinversion.  Finally, 
the  spline  coefficients  are  computed  via 

£  -  (L-1)1^.  (118) 


The  end  result  of  this  procedure  is  a  spline  solution  having  the  property  of 
"minimum  curvature." 


Several  pertinent  remarks  should  he  made  in  closing.  First,  there  is  no 
need  to  have  the  entire  matrix  [A:b]  in  computer  storage  at  one  time.  A 
sequential  Householder  procedure  [Lawson  and  Hanson  (1974)]  can  be  employed 
to  convert  [A:b]  to  an  upper  triangular  form;  only  m^  rows  of  data  need  to  be 
operated  on  at  a  time,  requiring  at  most  +  N  +  1  rows  of  storage.  The 
Integers  m^  can  be  as  small  as  1,  which  permits  the  greatest  economy  of  storage. 
The  resulting  least-squares  problem 


(119) 


PREVIOUS  PACE 
IS  BLANK 


where  R  is  upper  triangular,  is  equivalent  to  (115)  in  that  it  has  the  same 

set  of  solutions  and  the  same  minimal  residual  norm.  Furthermore,  the  matrix 

R  has  the  same  set  of  singular  values  as  the  matrix  A.  The  N  x  N  reference 

-1  T 

frame  t r  nsformation  (L  )  can  then  be  applied  to  the  R  matrix  to  yield 

RCL-1)1^  =  "d  (120) 

and  the  SVD  can  be  applied  to  this  result.  Since  the  matrix  L  *  is  lower 
triangular,  it  can  be  read  into  the  lower  triangular  portion  of  R  (plus  one 
additional  row  of  storage)  and  a  packed  multiply  can  be  performed. 

Second,  the  curvature  matrix  C  is  easily  computed  using  the  weighting 
matrices  W(n),  Y(n),  Z(n)  and/or  WP(n),  YP(n)  ZP(n). 

Third,  since  L  is  lower  triangular,  L  *  can  be  computed  strictly  through 
back  substitutions. 

Fourth,  the  spline  fit  is  controlled  through  specification  of  the  nodal 
spacing  and  the  signal-to-noise  ratio.  A  good  rule  of  thumb  is  to  have  the 
nodal  spacing  in  each  direction  equal  to  one-fourth  of  the  shortest  wavelength 
present  in  that  direction.  This,  in  part,  must  be  done  through  appropriate 
scaling  of  the  independent  variables  of  the  data. 

Finally,  if  an  M-variate,  minimum-curvature  spline  is  applied  to  M+l 
randomly  spaced  data  points,  an  M-dimensional  hyperplane  is  obtained.  While 
this  last  observation  has  little  practical  application,  it  does  nicely  illus¬ 
trate  the  minimum-curvature  property  of  these  multivariate  splines. 
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Figure  12.  Subset  of  Bivariate  Basis 
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APPENDIX  A 

The  Weighting  Matrices  W(n) 


defined  via 


The  weighting  matrices  W(n) ,  n-2,  3,  .  .  . 


APPENDIX  C 

The  Weighting  Matrices  Y(n) 


The  weighting  matrices  Y(n),  n«2,  3,  .  .  .  defined  via 

/  v  i= 1 »  • 


Y(n)  - 


m*j  <t>l  (t)dt 


are  given  below. 
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